- standard linear algebra algorithms for double and
complex values in Xi-
Define an example vector and matrix and transpose the matrix
( 1)>A={ {1,3}, {2,4} }; B=transpose(A); ( 2)>b={1,5};To perform a matrix multiplication use the "#" operator
( 3)>print(A#A); <intarr> 7 15 10 22The linear equation A*x=b will be solved by
( 4)>print(linear_solve(A,b)); <dblarr> 5.5 -1.5The command
( 5)>print(invert(A)); <dblarr> -2 1.5 1 -0.5gives the inverse of the matrix A. Furthermore to compute the determinant and to estimate the reciprocal condition numbers (first to the 1-norm, second to the infinity-norm) type
( 6)>print(determ(A),rcond(A),rcond_inf(A)); <double> -2 <double> 0.047619 <double> 0.047619Just as a reminder: rcond(A)=1 / ( norm(A) * norm(inv(A)) ). The eigenvalues and eigenvector can be obtained by
( 7)>print(eigenvalue(A),eigenvector(A)); <cpxarr> ( 5.3722813, 0) (-0.37228132, 0) <dblarr> 0.56576746 -0.90937671 0.82456484 0.41597356The eigenvectors are stored in the columns of the result matrix. The command
( 8)>[LU, pivot]=lu(A); ( 9)>print(LU, pivot); <dblarr> 2 4 0.5 1 <intarr> 2 2performs a LU decomposition of the matrix A. This is the first example where a function has more than one (here two) return arguments. The first return value LU contains the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored. The second one represents the pivot indices; row i of the matrix A was interchanged with row pivot[i]. Proceed this computation you can evalute the inverse of A by
( 10)>inv=lu_invert(LU, pivot);If you do not need the pivot indices, call the function by
( 11)>LU=lu(A);A QR decomposition can be computed by,
( 12)>[t, R]=qr(A); <dblarr> 1.4472136 0 <dblarr> -2.236068 -4.9193496 0.61803399 -0.89442719where the first return array has the scalar factors of the elementary reflectors as entries and the second value contains the Q and R matrix (for further details see the reference guide). Another supported decomposition is the SVD.
( 13)>[w,u,vt]=svd(A);Here w is a two element vector of singular values, u a (2 by 2) matrix and vt a (2 by 2) orthogonal matrix. A can be rewritten as
( 14)>W={ {w[0],0} , {0,w[1]} }; ( 15)>print(u#W#vt); /* u#W#vt equal to A */ <dblarr> 1 3 2 4Uses "back substitution" to solve a linear equation A*x=b
( 16)>x=svd_back_sub(w,u,vt,b);A typical use of the functions svd and svd_back_sub is to zero the very small w[i]'s before executing the back substitution to make the result stable against small changes in b. (for further details have a look in a standard book on matrix numerics).